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ABSTRACT 

A new method of polarimetric calibration is presented in which the instrumental response is derived 
from regular observations of PSR J0437— 4715 based on the assumption that the mean polarized emission 
from this millisecond pulsar remains constant over time. The technique is applicable to any experiment 
in which high-fidelity polarimetry is required over long time scales; it is demonstrated by calibrating 
7.2 years of high-precision timing observations of PSR J1022+1001 made at the Parkes Observatory. 
Application of the new technique followed by arrival time estimation using matrix template matching 
yields post-fit residuals with an uncertainty-weighted standard deviation of 880 ns, two times smaller 
than that of arrival time residuals obtained via conventional methods of calibration and arrival time 
estimation. The precision achieved by this experiment yields the first significant measurements of the 
secular variation of the projected semi-major axis, the precession of periastron, and the Shapiro delay; 
it also places PSR J1022+1001 among the ten best pulsars regularly observed as part of the Parkes 
Pulsar Timing Array (PPTA) project. It is shown that the timing accuracy of a large fraction of the 
pulsars in the PPTA is currently limited by the systematic timing error due to instrumental polarization 
artifacts. More importantly, long-term variations of systematic error are correlated between different 
pulsars, which adversely affects the primary objectives of any pulsar timing array experiment. These 
limitations may be overcome by adopting the techniques presented in this work, which relax the demand 
for instrumental polarization purity and thereby have the potential to reduce the development cost of 
next-generation telescopes such as the Square Kilometre Array. 

Subject headings: methods: data analysis — techniques: polarimetric — polarization — pulsars: 
general — pulsars: individual (PSR J1022-I-1001, PSR J0437-4715) 



1. INTRODUCTION 

High-precision pulsar timing exploits the exceptional 
rotational stability of millisecond pulsars to measure or- 
bital dynamics in binary systems and perform unique 
tests of gen eral relativity in the strong field regime (e.g. 
Stairs 2006). The basic measurement in timing analyses 
is the pulse time-of-arrival (TOA), the epoch at which a 
fiducial phase of the pulsar's periodic signal is received 
at the observatory. The difference between the mea- 
sured TOA and the arrival time predicted by a best-fit 
model is known as the timing residual. Timing residu- 
als with a standard deviation of the order of 100 ns have 
been ach ieved for a growing numbe r of millisecond pul- 
sars (e.g. lLommenll200ll ; iHotan 2007), which has renewed 
both theoretical and experimental interest in the detec- 
tion of gravit ational radiation using a pulsar timing ar 



igation (|Kocz et al.M2010t iNita fc Garvll2010D . compensa- 
tion for the eff ects of propagation through the interstellar 
medi u m (e.g., You et al.T2007at iHemberger fc Stinebrind 
120081: IColes et al.l I2010D . statistical analysis of the 
stochastic nature of the pul s ar signal ([Manchester et al.l 
19751: ICordes fc Downsl ll985l; iRathnas ree&Rankin 1995: 



Oslowski et al.l 120111). and the characterization of tim ing 



ng . 

ray (PTA; e.g.lFoster fc Backerlll990tlJaffe fc Backer[2 003: 
Wvithe fc Loebll2003l:lJenet et al.l2006l:ISesana et al.l200S : 



Yardlev et al.ll201lUvan Haasteren et alj|201lir 



The sensitivity of a PTA and the confidence limits 
that may be placed on any gravitational wave detec- 
tion directly depend upon the precision and accuracy 
with which pulse arrival times can be estimated. There- 
fore, an important part of the ambitious PTA effort 
is the quantification and reduction of systematic error 
through the development of improved methods of ar- 
rival time estimation (e.g.. IHotan et al" 20051: Ivan Stratenl 
2006 ), instrumen tal calibration ( Jenet fc Anderson! 119981 : 



van Stratenll2004l ). radio-frequency interference (RFI) mit 



noise (ICordes fc Downslll985l: IShannon fc Cordesll2010D . 

At the Parkes Observatory, the majority of pulsar timing 
data are obtained from observations of PSR J0437— 4715, 
the closest and brightest milliseco nd pulsar known. Dis - 
covered in the Parkes 70-cm survey ([Johnston et al.llT993l) . 
PSR J0437-4715 has a spin period of ~ 5.7 ms; at 20 cm, 
it has a pulse width at half maximum of about 130 /is 
(Na varro et al.l 19971) and an average flux of 140 mJy 
(|Kramer et al.|[l998l ). Owing to its sharp peak and large 
flux density, it is an excellent target for high-precision pul- 
sar timing studies. It is also fortuitously located away 
from the Galactic plane (/ = 245°, b = -42°), which per- 
mits long observing sessions during periods of local sidereal 
time when competition for the telescope is relatively low. 

Early timing of this pulsar highlig hted the need for mor e 
accurate polarimetric calibration (|Sandhu et al.l 119971 ). 
The systematic timing error due to instrumental polariza- 
tion artifacts is readily observed in the dramatic variation 
of arrival time residuals as a function of parallactic angle 
(see Figure [IJ . This issue motivated the discov ery and 
devel opment of the polarimetric invariant profile ([Brittonl 
2000). Arrival times derived from the invariant profile 
proved to be significantly more accurate than those of the 
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Fig. 1. — Systematic timing error in one day of uncalibrated PSR J0437— 4715 observations. As a function of time, the receiver rotates with 
respect to the sky, thereby geometrically altering the instrumental response to the highly polarized pulsar signal. The induced systematic 
error is about an order of magnitude larger than the estimated arrival time uncertainty (~ 250 ns). 



total intensity profile, a breakthrough that led to the hrst 
detection of an nual-orbital parallax an d a new test of gen- 
eral relativity (jvan Straten et aT1l200lD . 

Despite its demonstrated advantages over the total in- 
tensity profile, there remain a number of drawbacks to 
the use of the invariant profile. First, the signal-to-noise 
ratio S/N of the invariant profile is lower than that of 
the total intensity, reaching zero for completely polarized 
sources. Second, the noise in the invariant profile is nei- 
ther normally distributed nor homoscedastic with respect 
to pulse longitude (see Appendix |A"|) : these two proper- 
ties are assumed, either explicitly or implicitly, by most 
template matching algorithms commonly used for pulsar 
timing. Finally, the computation of the invariant profile 
requires estimation and correction of the bias due to noise, 
and the relative precision of the bias estimate is inversely 
proportional to the square root of the number of discrete 
samples within off-pulse regions of pulse longitude. For a 
typical observation of PSR J0437— 4715, the relative error 
in the bias correction is of the order of 8%, an unaccept- 
ably large uncertainty when aiming for arrival time uncer- 
tainty that is four orders of magnitude smaller than the 
pulse width. 

These concerns motivat ed the developme nt of matrix 
template matching fMTM: lvan Stratenll2006| ). a technique 
that exploits the additional timing information in the po- 
larized component of the pulsar signal while simultane- 
ously eliminating residual polarimetric calibration errors. 
For some pulsars, the mean polarization profile exhibits 
sharper features than observed in the total intensity; these 
features generate more power at higher harmonics of the 
pulsar spin frequency, yielding greater arrival time preci- 
sion than might be exp ected from the pol arized flux alone 
(e.g. see Figure 1 of Ivan Stra ten 2006) . Furthermore, 
MTM improves arrival time accuracy by mitigating sys- 
tematic timing errors due to instrumental polarization ar- 
tifacts. 

Although the development of MTM was primarily 
driven by the hig h-precision timing of PSR J0437—4715, 
Ivan Straten! (|2006j) predicted that PSR J1022+1001 would 
benefit most from this technique, including a ~ 33% 



and a significant decrease 
PSRJ1022+1001 was con- 
temporaneously discovered in both the P r inceto n- Arecibo 
Declination-Strip Survey ()Camilo et aL 



increase in timing precision 
in systematic timing error. 



1996) a nd th e 



Green Bank North ern Sky Survey (|Saver et al.l Il997l ). 
iCamilo et al.l (|1996|) observed variations in the shape of 
its mean pulse profile on a timescale of the order of min- 
utes and noted that the accuracy of arrival time esti- 
mates was limited by this effect. Mean profile varia- 
tions over such timescales would be distinct from pulse-to- 
pulse shape variations such as giant pulses, microstructure, 
and d rifting s ub-pu lses (e.g. lHankinsI 119711 : iLvne et ail 
119711 : iBackerl I1973D . which when ave raged over many 
pulse s can be described as phase jitter (jCordes fc Downsl 
1985) or stoch astic wide-band imp ulse-modulated self- 
noise (SWIMS; lOslowski et~al"1l201ll) . As it is typically 
assumed that the mean pulse profile does not change with 
time, detection of systematic profile variations in a mil- 
lisecond pulsar would have a serious impact on the field of 
high-precision timing. 

Accordingly, the stability of the PSR J1022+1001 pulse 
profi le has been the sub j ect of a number of detailed stud- 
ies dKramer et all 119991 : iRamachandran fc Krame"rl [2003; 
iHotan et al.ll2004al) . These are reviewed in ^ following a 
mathematical treatment of the systematic timing error due 
to instrumental distortion of the total intensity profile. 
describes the observing system and the novel procedures 
used to calibrate the PSR J1022+1001 data for this study, 
including a new scattered power correction algorithm that 
is described in more detail in Appendix [Bj As part of this 
analysis, a 7.2-year history of the polarimetric response at 
Parkes is presented and its impact on systematic timing 
error is discussed. In 21 the results of this study are com- 
pared and contrasted with previous work, the relevance of 
the calibration method to other pulsars is discussed, and 
potential future directions are considered. 

2. SYSTEMATIC TIMING ERROR 

In the following discussion, the polarization of electro- 
magnetic radiation is described by the second-order statis- 
tics of the transverse electric field vector e, as represented 
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by the four Stokes parameters 

S fc = <eW>- (1) 
The angular brackets denote an ensemble average; e* is 
the conjugate transpose of e; <x is the 2x2 identity ma- 
trix; cr i, er 2 , and er 3 are the Pauli spin matrices; So is 
the total in tensity, and S = (Si, S 2 , S3) is the polarization 
vector (e.g. lBrittorj|2000f) . The reception and propagation 
of polarized radiation is described by linear transforma- 
tions of e, as represented using complex 2x2 Jones ma- 
trices. Any non-singular matrix can be decomposed into 
the product of a unitary matrix and a self-adjoint matrix. 
Unitary transformations result in a three-dimensional Eu- 
clidean rotation of the Stokes polarization vector, leaving 
the total intensity unchanged and preserving the degree of 
polarization. Self-adjoint Jones matrices effect a Lorentz 
boost of the Sto kes 4- vecto r, mixing total intensity and 
polarized flux (|Brittonll2000|). The boost c omponent, also 
known as the voldistortion (jHa makcr 2000h . of the instru- 
mental response distorts the shape of the total intensity 
pulse profile and introduces systematic arrival time errors. 
Physically, instrumental boost transformations are caused 
by differential amplification of the signals from the two 
orthogonally-polarized receptors and by cross-coupling be- 
twee n these recep t ors. 

In Ivan Stratenl ([2006) , the systematic timing error in- 
duced by a given level of polarization distortion is pre- 
dicted for a number of millisecond pulsars using a numeri- 
cal simulation in which the template profile for each pulsar 
is copied and subjected to a boost transformation. The 
orientation of the boost axis in Poincare space is varied, 
and the minimum and maximum induced phase shifts be- 
tween the distorted total intensity profile and the template 
profile are recorded. It is also possible to predict the sys- 
tematic timing error due to polarization calibration errors 
by computing the rate at which the best estimate of the 
phase shift varies with respect to the instrumental boost 
tr ansformation param eters. Beginning with equation (12) 
of Ivan Stratenl ([2006D . the best estimate of the phase shift 
ip derived using only the total intensity is given by 

I So 

Here, So,m = S'$ (v m )So(v m ) is the cross-spectral power 
of the observed total intensity S' and the template total 
intensity So in the Fourier domain, <f>Q m is the argument 
of So mi and the summations in the numerator and de- 
nominator are performed over m = 1 to m = N, where N 
is the maximum harmonic at which the fluctuation power 
spectrum exhibits significant power. 

Now consider an observed total intensity profile that 
is a copy of the template subjected to an instrumental 
boost transformat ion, as parameterized in Section 4.1 of 
Ivan Stratenl (|2004h such that sinh 2 /3 = b • b. Referring 
to equation (10) of lBritto"nl (pOO O), the transformed total 
intensity is given by 

S = (1 + 2b-b)S + 2b a b ■ S, (3) 
where bo — cosh/3 = \/l + b • b. The partial derivatives 
of S' with respect to the boost parameters bk are used 
to compute the partial derivatives of the phase shift with 
respect to these parameters when the boost parameter /3 
is zero, 



•Pk 



dip 
db k 



/9=0 



7r£ |S (^m)| 2 z4 



The above expression defines the components of a three- 
dimensional gradient ip. The magnitude of the gradient 
•pp = \<p\ provides a measure of the susceptibility of ar- 
rival time estimates from a given pulsar to instrumental 
distortion. 

To aid in the physical interpretation of the susceptibil- 
ity tpp, the first order approx imations to equations (25) 
and (32) of Ivan Stratenl (2006) are used to define the sys- 
tematic timing error induced by either a 1% differential 
gain error or receptor non-orthogonality of 0.01 radians 
(~ 0.6°), 



t p = 5 x 10-V^P, 



(5) 



where P is the pulsar spin period. Values of for 
each of the 20 millisecond pulsars that are regularly ob- 
served as part of the Pa rkes Pulsar Timing Array (PPTA; 
iManchester et al.l 12012ft project are reported in Table [TJ 
The quantities in this table are der ived from the hig h S/N 
polarization profiles presented by lYan et al.l ([201 If ). For 
the majority of pulsars in the PPTA, a modest degree of 
instrumental distortion induces systematic timing error of 
the order of 100 ns, which is sufficient to seriously impede 
progress toward the det ection of the sto chastic gravita- 
tional wave background ([Jenet et al.ll2005l ). 

To first order, the systematic timing error due to in- 
strumental distortion is given by the inner product, b • (p. 
In the polar coordinate system best suited to describing 
linearly polarized receptors, b ~ (7, dg,§ e ), where 7 pa- 
rameterizes the differential receptor gain, and 5 e and Sg 
describe t he non-orthogonal ity of the receptors (see Sec- 
tion 3.3 of Ivan Stratenl [20061 ). These properties of the in- 
strument may vary as a function of time for a variety of 
reasons. Diurnal variations are introduced by the paral- 
lactic rotation of the receiver with respect to the sky, as 
shown in Figure [TJ Furthermore, step changes in differ- 
ential gain are introduced when the levels of attenuation 
applied to each of the two orthogonal polarizations are 
reset at the start of each observation and/or periodically 
updated during the observation. In contrast, the cross- 
coupling of the receptors is typically assumed to remain 
relatively stable over timescales of the order of months. 

Of particular interest to a long-term timing experiment 
such as a pulsar timing array are any step changes in b due 
to modifications of the instrumentation and/or slow vari- 
ations in b due to gradual degradation of one or more sys- 
tem components. The temporal variations in distortion- 
induced error due to the long-term evolution of a given 
instrument will be correlated between the different pul- 
sars observed using that system. In an array of N pulsars, 
the spectral structure of the correlated systematic error 
due to instrumental distortion is described by 



C T (/) = TCp(/)T a 



(6) 



where C T is the NxN matrix of the auto- and cross-powe r 
spectra of the timing residuals (e.g. lYardlev et al.|[20lT[ ). 
Y is the N x 3 Jacobian matrix, 



dTj_ 
db k 



P 



dip 



3 db 



(7) 



(4) 



(Pj is the spin period of the jth pulsar) and Cp(f) is the 
3x3 matrix of the auto- and cross-power spectra of the 
instrumental distortion parameters (the components of b). 
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Table 1 

Relative Arrival Time Uncertainties and Systematic Timing Error 



Pulsar 


a 


■■p 


a 


$ 


Tfj (ns) 


J0437-4715 





,85 


1 


.43 


207 


J0613-0200 





,92 


1 


.46 


59 


J0711-6830 





,88 


1 


.54 


81 


J1022+1001 





.68 


1 


.65 


282 


J1024-0719 





.74 


2 


.11 


34 


J1045-4509 





.88 


1 


.48 


338 


J1600-3053 





.90 


1 


.39 


115 


J1603-7202 





.85 


1 


.55 


142 


J1643-1224 





.91 


1 


.40 


266 


J1713+0747 





.85 


1 


.58 


6 


J1730-2304 





.71 


1 


.70 


198 


J1732-5049 





,96 


1 


.38 


185 


J1744-1134 


1 


.56 


6 


.43 


105 


J1824-2452A 





,88 


2 


.56 


18 


J1857+0943 





.89 


1 


.43 


124 


J1909-3744 


1 


.02 


1 


.51 


22 


J1939+2134 





.95 


1 


.49 


44 


J2124-3358 





.85 


1 


.45 


127 


J2129-5721 


1 


.15 


1 


.61 


211 


J2145-0750 





.95 


1 


.44 


147 



In practice, an optimally-weighted sum of the cross- 
power spectral harmonics is used to define a detec- 
tion statistic with maximum sensitivity to a stochas- 
tic background of low-frequency gravita tional waves (e.g. 
lAnholm et all 120091: lYardlev et all 1201 ll) . Given only the 
mean polarization profile of each pulsar in the timing ar- 
ray, it is possible to predict the impact of correlated sys- 
tematic timing error due to instrumental distortion on 
such a detection statistic. First, only the dominant har- 
monic of the cross-power spectrum of two pulsars with 
distortion gradients ip A and tp B is considered. Next, it 
is assumed that the instrumental distortion parameters 
are uncorrelated and homoscedastic, such that C T — cr^I, 
where I is the 3x3 identity matrix. Following these sim- 
plifications, a first-order approximation of the coefficient 
of correlated systematic timing error is given by 

cab = | . || . | - (8) 

This correlation coefficient depends only on the shape of 
the mean polarization profile and has no modal struc- 
ture on the celestial sphere. That is, the correlated tim- 
ing error due to instrumental distortion will corrupt all 
of the moments in a multipole expansion of timing de- 
lay, which will impact on all of t he primary goals of an y 
pulsar timing array experiment (|Foster fc Backer! [T990) . 
including the long-term measur ement of time (monopole 
moment; iGuinot &: Petit! |1991|) . correct ions to the So- 
lar S ystem ephemeris (dipole moment; IChampion et all 
|2010| ). and detection o f the gravitational wave background 
(quadrupole moment; iHellings fc Downslll983|) . 

To illustrate the potential importance of correlated sys- 
tematic timing e rror, the h igh S /N polarization pro- 
files presented by lYan et all (|2011l) are used to compute 
the cross-correlation coefficients for each of the 15 pul- 



sar pairs presented in Figure 4 of lYardlev etail (|20TTI). 
The c ross-spectral power measurements of lYardlev et al.l 
(|201l[) are compared with the predicted values of cab 
in Figure the coefficient of correlation between pre- 
dicted and measured values is p ~ 0.35. The agreement 
between predicted and experimental measures of corre- 
lation indicates that the systematic error due to instru- 
mental distortion is a plausible contributor to the appar- 
ent an ti-correlation between the results of lYardlev et al.l 
(1201 ll) and the expec t ed qu adrupolar signature predicted 
bv IHellings fc Downs! (|1983l) . 

Also shown in Table Q] are the relative arrival tim e un- 
certainties a v and defined in Ivan Stratenl (120061. The 
former is the ratio between the theoretical error in arrival 
time estimates derived using matrix template matching 
(MTM) and the error in estimates derived using only the 
total intensity; the latter is the ratio between the uncer- 
tainties in arrival times derived from the invariant interval 
and from the total intensity. The statistical uncertainty 
in the invariant profile arrival times is underestimated be- 
cause the effects of heteroscedasticity and bias correction 
error are not considered. Regardless, for every pulsar in 
the PPTA, precision is lost through use of the invariant 
interval. The increase in uncertainty is most dramatic for 
PSR J1744-1134, which emits almost 100% linearly polar- 
ized radiation. For this pulsar, MTM is also predicted to 
yield arrival times with greater error than those derived 
from the total intensity owing to the large coefficient of 
multiple correlation (~ 0.9) between the phase shift used 
to compute the TOA and the model parameters that de- 
scribe the polarization transformation. 
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Fig. 2. — Comparison between predicted and measured estimates of correlation in pulsar timing residuals. The predicted coefficients of 
correlated systematic timing error, cab (cross es) are qualitat i vely similar to the optimally-weighted sum of cross-spectral power estimates 



A?-£(0y) (points with error bars) presented bv lYardlev et all 1)20111 ). The values of cab have been scaled by 2.7 X 10 , as determined by 
a (non- weighted) least-squares fit to the values of Ay£(0ij). 



2.1. PSRJ1022+1001 

Table Q] also predicts that matrix template matching 
will yield the greatest relative improvement in arrival time 
precision for PSR J1022+1001. This pulsar also has the 
second-highest value for predicted systematic timing er- 
ror owing to the susceptibility of the total intensity profile 
to polarization distortion. This may explain the profile 
sh ape variations and e xcess arrival time error first noted 
by ICamilo et al.l (|l996t ). In a detailed study of the spec- 
trum and polarization of thes e variations using the Effels- 
berg 100-m Radio Telescope, iKramer et al.l (|1999( ) assert 
that the observed fluctuations in pulse shape cannot be 
explained by instrumental effects. However, the data used 
in this analysis were calibrated using the "ideal feed" as- 
sumption that there is no cross-coupling between recep- 
tors and that the reference source used for gain calibra- 
tion produces an equal and in-phase response in each re- 
ceptor. Given that receptor cross-coupling of the order 
of 10% to 20% is commonly obse rved when more accu- 
rate calibration is performed (e. g. IStinebring et al.l 11984 
IXilourisI Il99lb Ivan Stratenl 12004 ) , an instrumental origin 
of the variability deserves consideration. 

For instance, the temp oral variation of the total in- 
tensity profile observed by IKramer et al.l Jl999) is repro- 
ducible using a simple model in which the polarized signal 
is rotated by the parallactic angle of the receiver and then 
subjected to an instrumental boost that mixes ~10% of 
the linearly-polarized flux with the total intensity. The 
similarity between the simulated observations plotted in 
Figure Bland the total intens ity profiles presented in Fig- 
ure 2 of Kramer et al.l (|1999[ ) supports the argument that 
the variations are the result of uncalibrated instrumental 
distortion. 

Similarly, the small f requency-scale variations observed 
by IKramer et al.l (|1999| ) may be readily explained by spec- 
tr al variation of the instrumental response (e.g. Figure 3 
of Ivan Stratenl 12004 . The observations used to demon- 
strate th e narrow-band variat ions and presented in Fig- 
ure 8 of IKramer et al.l (fl999 ) were recorded using the 



Effelsberg-Berkeley-Pulsar-Processor with a bandwidth of 
56 MHz. In this mode, only the fluxes in each of the 
circularly polarized signals (i.e. two out of four Stokes pa- 
rameters) are retained; consequently, only the mixing be- 
tween circularly-polarized flux and total intensity can be 
calibrated in these data. Therefore, frequency-dependent 
mixing between linearly-polarized flux and total intensity 
is a plausible explanation for the observed profile shape 
v ariations. 

iRamachandran fc KrameU (|2003f ) analyzed the variabil- 
ity of PSR ,11022+1001 using multi-frequency observa- 
tions made at the Westerbork Synthesis Radio Telescope 
(WSRT). They argue that, because the WSRT is com- 
posed of equatorially-mounted antennas, it is possible to 
rule out pulse shape variations due to instrumental dis- 
tortion combined with parallactic rotation of the feed. 
However, the polarimetric response of a phased array can 
vary in both time and frequency for a wide variety of rea- 
sons, including imperfect phase calibration and (as for any 
single-dish antenna) instrumental gain variations. Indeed, 
the automatic gain controllers us ed to prevent instrumen- 
tal saturation (|Voute et al.l 12002]) also make it impossible 
to accurately calibra te the polarimetric respon se of the 
WSRT phased array (|Edwards fc Stapperdl200l . 

The high degree of susceptibility of the PSR J1022+1001 
pulse profile to polarization distortion and the inadequate 
instrumental calibration employed in these previous works 
cast doubt on the conclusion that the profile variat i ons are 
intrinsic to the pulsar. Furthermore, iHotan et al . (2004a) 
find no evidence of significant profile shape variations in 
PSR J1022+1001 observations made at the Parkes Obser- 
vatory. This study achieved arrival time residuals with a 
standard deviation of the o rder of 2.3 us and y 2 per degree 
of freedom ~ 1.4, whereas IKramer et al.l ( 1999h obtained 
a reduced x 2 close to unity only after systematic errors of 
the order of 10 /is were added in quadrature to the formal 
uncertainties. 

The analysis bv lHotan et al . (2004a) is based on obser- 
vations integrated over 5 min and 64 MHz and spanning 
1.3 years. Increasing the integration length to 60 min and 
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Fig. 3. — Simulated temporal variation of the PSR J1022+1001 mean total intensity profile. The parallactic an gle of each simulate d 
observation is computed for the Effelsberg 100-m Radio Telescope at the epoch of the corresponding panel in Figure 2 of lkramer et al.l l|1999f l. 
Over this 3.3 hr observation, the parallactic angle varies by 37°. 
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the time span to 2.3 years yielded arrival time residuals 
with a standa rd deviation of the order of 1.5 /is and re- 
duced x 2 ~ 3 (jHotan et al.ll2006t) . With t he same instru- 
menta tion and similar integration lengths, IVerbiest et al.l 
(2009) achieved 1.6 /is residuals in data spanning 5.1 years. 
In the absence of systematic error and other sources of 
noise (e.g. dispersion variations), the one- hour integra- 
tions used in these two studies should have yielded timing 
residuals of the order of 660 ns. This is of the same order 
as the systematic timing error introduced by a differential 
gain error of ~ 2.3% or by receptor non-orthogonality of ~ 
1.4° (cf. TablejT]). However, the arrival times used in these 
studi es were derived from either u ncalibrated dHotan et ail 
2006| or inaccurately calibrated (jVerbiest et al.ll2009D to- 
tal intensity pulse profiles, which are highly susceptible 
to instrumental distortion. Therefore, it is reasonable to 
anticipate that high-fidelity polarimctric calibration will 
improve the accuracy of PSR J1022+1001 arrival time es- 
timates. 

This conjecture is verified in the following section, which 
describes a new method of deriving the instrumental re- 
sponse by using PSR J0437— 4715 as a polarized reference 
source. The method is used to calibrate high-precision 
timing observations of PSR J1022+1001 and demonstrated 
to reduce the systematic error due to instrumental polar- 
ization distortion. 

3. OBSERVATIONS AND ANALYSIS 

Dual-polarization observations of PSR J0437— 4715 and 
PSR J1022+1001 were made at the Parkes Observatory us- 
ing the H-OH receiver and th e center element of the Parke s 
21-cm Multibeam receiver (jStavelev-Smith et al.l Il996h . 
Two 64-MHz bands, centered at 1341 and 1405 MHz, were 
digitized using two bits per sample and processed using 
the second genera tion of the Caltech-Par kes-Swinburne 
Recorder (CPSR2: [Baiiesl[2003t lHotan|[2007h . To maintain 
optimal linear response during digitization, the detected 
power was monitored and the sampling thresholds were 
updated approximately every 30 se conds. The baseband 
data were reduced using DSPSR (jvan Straten fc Bailesl 
1201 If) , which corrects quantization distortions to the 
voltag e waveform using the d ynamic level setting tech- 
nique (jJenet fc Andersonll 19981 hereafter J A98), then per- 
forms phase-coherent di spersion removal (jHankins! 1 19711 : 
Han kins fe Rickettl 1975) whi le synthesizing a 128-channel 
filterbank (jvan Stratenll2003l) . The Stokes parameters are 
then formed and integrated as a function of topocentric 
pulse phase, producing uncalibrated polarization profiles 
with 1024 discrete pulse longitudes. 

3.1. Minimizing quantization distortions 

During analog-to-digital conversion, the radio signal is 
subjected to non-linear distortion s that significantly alter 
the o bserved pulse profile (JA98; iKouwenhoven fc Voutd 
2001). JA98 introduce a dynamic output level setting 
technique that is employed by DSPSR to correct the un- 
derestimation of undigitized power. To implement this 
technique, the digitized data are divided into consecutive 
segments of L samples and, for each segment, the number 
of low-voltage states M is counted. A histogram of oc- 
currences of M is archived with the pulsar data and used 

1 http : //psrchive . sourcef orge . net 



offline to evaluate the quality of the digital recording. 

When the voltage input to the digitizer is normally dis- 
tributed, the ratio $ = Al/L has a binomial distribution 
as in equation (A6) of JA98. The difference between this 
theoretical expectation and the recorded histogram of M 
provides a measure of the degree to which either the input 
signal is not well described by a normal distribution or the 
sampling thresholds diverge from optimality (e.g. at the 
start of an observation). This difference, called the two-bit 
distortion, is given by 



D=J2 i v ( M / L ) - H{M)] 



(9) 



where V{$) is the expected binomial distribution and 
H(M) is the recorded distribution of M. Separate his- 
tograms of M are maintained for each polarisation, and 
the reported distortion is simply the sum of the distortion 
in each polarisation. 

For each 16.8-second integration output by CPSR2, the 
two-bit distortion D was computed and all data with 
D > 3.5 x 10~ 4 were discarded. This threshold was chosen 
using the histograms plotted in Figure 01 beyond the range 
of plotted values, a sparse tail of outliers extends to a maxi- 
mum value of D ~ 0.98. Approximately 3.3% of the obser- 
vations were rejected; the remaining data were corrected 
for scattered power using the method described in Ap- 
pendix[B]and implemented in the psrchive s oftware pack- 
age fo r pulsar data archival and analysifQ (jHotan et al.l 
2004b|). The corrected data were then summed until the 
integration length was of the order of five minutes. 

3.2. Polarimetric template profile 

To use a pulsar as a polarized reference source requires 
a high S/N, well-calibrated, template polarization profile. 
From over 700 hours of PSR J0437— 4715 observations were 
selected only those sessions with low levels of two-bit dis- 
tortion and of sufficient duration to achieve a precise es- 
timate of the instrumental response. About 50 sessions, 
each 8 hours in duration and with corresponding pulsed 
noise diode observations, were c alibrated using me asure- 
ment equation modeling (MEM; Ivan Stratenll200l . The 
observed Stokes parameters were normalized to account 
fo r pulsar flux variat ions as described in Appendix As 
in Ivan Straten! (|2004[) , the two degenerate model parame- 
ters (a boost along the Stokes V axis and a rotation about 
this axis) were constrained by assuming that observations 
of 3C 218 (Hydra A) have negligible circular polarization 
and that the orientation of the receiver is known. Although 
these assumptions are not necessarily valid, absolute cer- 
tainty in the Stokes parameters is not required for the 
purposes of a high-precision timing experiment. It is suffi- 
cient that the observed Stokes parameters can be mapped 
to the template profile by a single Jones transformation. 

This prerequisite may be evaluated using the objective 
merit function of the MTM least-squares fit; the reduced 
Y 2 statistic will be greater than unity whenever the trans- 
formation between the template and the observed pulse 
profile is not well-described by a single Jones matrix. Var- 
ious phenomena may contribute to a poor model fit. For 
example, if the mean polarized emission from the pulsar 
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Fig. 4. — Stacked histograms of the two-bit distortion D in the CPSR2 bands centered at 1341 MHz (top/gray) and 1405 MHz (bottom/light 
gray). 



evolves over time, then the difference between the tem- 
plate and the observed profile will vary as a function of 
pulse longitude in a manner that cannot be adequately 
modeled by a single Jones transformation. Furthermore, 
as Jones matrices describe only linear transformations of 
the electric field, matrix template matching will fail to 
model any non-linear component of the instrument. 

Analog-to-digital conversion using only two bits per 
sample is an intrinsically non-linear process; therefore, 
variations in the response of the digitizer adversely affect 
the MTM goodness-of-fit, as shown in Figure [5] Here, 
5 of the - 50 MEM-calibrated PSR J0437-4715 observa- 
tions were integrated to produce a template profile and the 
mean polarization profiles from each of the original ~ 50 
observations were fit to this template using MTM. The 
reduced x 2 °f eacn MTM fit (averaged over all frequency 
channels) is highly correlated with the median value of the 
two-bit distortion D in each session; i.e. greater distortion 
reduces the merit of the MTM fits. For the 5 observations 
used to create the template profile, x 2 /Nf Tee < 1 because 
covariance between the observation and the template is 
not included in the definition of the MTM objective merit 
function. 

Figure [5] justifies the experimental design decision to 
regularly update the two-bit sampling thresholds during 
each observation. Independent sampling thresholds were 
applied to each of the two orthogonal polarizations; con- 
sequently, the astronomical signal was subjected to an un- 
known differential gain that varies with time. Such differ- 
ential gain variations can be modeled using MTM, whereas 
the non-linear response of the digitizer cannot. Therefore, 
when recording the baseband signal with a two-bit digi- 
tizer, it is more important to maintain optimal linear re- 
sponse than it is to maintain a constant flux scale. 

Out of the ~ 50 sessions that were calibrated using 
MEM, 16 were selected with low median two-bit distor- 
tion and MTM reduced x 2 within 4% of unity. The MTM- 
corrected mean profiles from these 16 sessions were inte- 
grated to form a polarimetric template profile with a total 
integration length of ^100 hours, shown in Figure |5] The 
over-polarization seen from pulse phase 0.88 to 0.94 in this 
figure is an instrumental artifact that is currently not un- 
derstood. It may be that the unpolarized quantization 



noise is overestimated by the two-bit scattered power cor- 
rection algorithm described in Appendix [Bj However, a 
very similar artifact is observed in polarization data ob- 
tained with the seco nd-generation of the Parkes Digital 
Filter Bank (PDFB2; lYan et alj|20lTl ). which employs an 
8-bit digitizer. These facts might point to the existence 
of a non-linear component in the signal path shared by 
CPSR2 and PDFB2. However, the third and fourth gen- 
erations of the PDFB also regularly exhibit ~ 5% over- 
polarization when it is not observed using the latest gen- 
eration of baseband recording instrumentation at Parkes. 
Therefore, it may be only coincidence that a similar arti- 
fact is observed in both PDFB2 and CPSR2 data. 

3.3. Measurement Equation Template Matching 

Matrix template matching has been previously utilized 
to calibrate and derive arrival tim e estimates from ob - 
servations of PSR J0437— 4715 (e.g. IVerbiest et al.ll2008D . 
However, this method cannot be applied directly to 
PSR J1022-I-1001 because the mean flux o f this pulsar at 
20 cm (~6 mJy; Manche ster et al] 120121 ) is much lower 
and the precision of calibrator solutions derived from sim- 
ilar integration lengths is inadequate. It is not possible 
to integrate over longer intervals because the instrumen- 
tal response varies with time, primarily due to the par- 
allactic rotation of the receiver. Therefore, to calibrate 
the PSR J1022+1001 data, a new technique was devel- 
oped that combines MTM with MEM. The new method, 
called Measured Equation Template Matching (METM), 
matches multiple pulsar observations to the template pro- 
file, thereby increasing the precision of the solution. As 
with MEM, variations of the instrumental response (e.g. 
the parallactic rotation of the receiver) are included in the 
model. METM also incorporates measurements of an ar- 
tificial reference source (e.g. a noise diode coupled to the 
receptors), enabling the backend comp onent of the cali- 
brator solution to be later updated as in lOrd et all (|2004ft 
so that the solution may be applied to other astronomi- 
cal sources. In contrast with MEM, the METM method 
does not require any additional observations of a source of 
known circular polarization to constrain the boost along 
the Stokes V axis because the template profile provides 
this information. Using the new METM technique, a de- 
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Fig. 5. — Average matrix template matching goodncss-of-fit \ 2 /^free as a function of median two-bit distortion D in the CPSR2 band 
centered at 1341 MHz. 
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Fig. 6. — Mean polarization of PSR J0437— 4715, plotted as a function of pulse phase using polar coordinates: cllipticity, e, orientation, S, 
and polarized intensity, S = \S\ (plotted in grey below the total intensity, So). Flux densities are normalized by <rrj, the standard deviation 
of the off-pulse total intensity phase bins. Data were integrated over a 64 MHz band centered at 1341 MHz for approximately 96 hours. 
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Fig. 7. — Temporal and spectral variation of receptor cross-coupling. The non-orthogonality parameters, 8g (upper panel) and S e (lower 
panel), describe the mixing between total intensity and polarized flux in the 21-cm Multibeam receiver and down-conversion system of one 
of the CPSR2 bands. These estimates were derived from regular timing observations of PSR J0437— 4715 using the new METM method 
presented in this paper. The two gaps in the data correspond to periods of maintenance on the 21-cm Multibeam receiver. 
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Fig. 8. — Temporal and spectral variation of the non-orthogonality parameters 5g (black) and <5 E (grey). From top to bottom, the panels 
plot the solutions obtained on 11, 12, 13, 14 and 16 February 2010. The length of each error bar is twice the formal standard deviation 
returned by the least-squares fit. 



tailed history of the instrumental response at Parkes was 
derived from the PSR J0437— 4715 data and used to cali- 
brate the PSR J1022+1001 data. 

First, the high S/N template profile (created as de- 
scribed in the previous section) and the new METM 
method were used to derive a calibrator solution from ev- 
ery observation of PSR J0437— 4715 with more than 3.5 
hours of integration. A total of ~ 350 solutions were 
produced, spanning 2003 April to 2010 June. The in- 
strumental res ponse was modeled using equation (19) of 
iBrittonl (|2000D . which includes three independent Lorentz 
boost transformations described by the differential gain 
and the non-orthogonality parameters, S e and 6$. The 
derived non-orthogonality parameter estimates are plot- 
ted in Figure [7] In addition to smooth variations in Sg 
and 5 e , there are a number of distinct steps that corre- 
spond to physical changes in the 21-cm Multibeam receiver 
and downconversion system at Parkes. The gaps in the 
data correspond to maintenance periods during which the 
Multibeam receiver was unavailable (2003 November 9 to 
2004 September 7, and 2006 December 16 to 2007 May 
17). Immediately following the first maintenance period, 
the receptor cross-coupling remains stable for a period of 
approximately one year. The distinct changes in state that 
occur around MJD 53660 (2005 October 17) and MJD 
53800 (2006 March 6) are due to changes in the down- 
conversion system (the signal paths for the two CPSR2 
bands were swapped). Following the second maintenance 
period, the cross-coupling parameters are notably less sta- 
ble in both time and frequency. Figure [8] provides a closer 
look at the variations of these parameters over a period of 
one week and a bandwidth of 20 MHz. In addition to sig- 



nificant quasi-periodic variations in these parameters as 
a function of radio frequency, the spectral structure ap- 
pears to drift as a function of time. For example, from 11 
to 12 February 2010, the main features in the spectrum 
have shifted toward lower frequencies by ~ 1.5 MHz. Af- 
ter MJD ~ 54300, it is invalid to assume that the receptor 
cross-coupling remains stable for any period of time longer 
than a day. 

Figure[7]plots the non-orthogonality parameters in units 
of centiradians; these values can be multiplied by the es- 
timates of Tp listed in Table Q] to approximate the vari- 
ation in systematic arrival time error introduced by re- 
ceptor cross-coupling. For example, for PSR J1022+1001, 
the peak-to-peak variations of ~ 0.06 radians in both Sg 
and S e are predicted to induce systematic timing errors of 
the order of 1.7 /its. However, after integrating over all fre- 
quency channels, the net effect of the cross-coupling vari- 
ations will be diminished. Furthermore, receptor cross- 
coupling describes only two of the three forms of instru- 
mental distortion. In addition to the presumption of zero 
receptor cross-coupling, calibration techniques based on 
the ideal feed assumption typically also incorporate the 
false premise that the artificial reference source illuminates 
both receptors identically. Under this assumption, any in- 
trinsic imbalance in the induced amplitude of the reference 
source signal in each receptor will be misinterpreted as dif- 
ferential gain. In contrast, the METM technique includes 
direct estimation of th e Stokes paramete rs of the reference 
source (as in Fig. 2 of Ivan Straten|[2004T) . Using these es- 
timates, the systematic error in differential gain 5 7 due to 
unbalanced reference source amplitudes may be derived. 

Given estimates of <L, Sg and S e as a function of fre- 
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Fig. 9. — Systematic arrival time error induced by instrumental distortion in the 21-cm Multibeam and H-OH receivers at Parkes. The 
periods during which the H-OH receiver was used are marked by two horizontal lines near the MJD axis. 



qucncy, it is possible to directly compute the systematic 
arrival time error at each epoch by applying the instru- 
mental distortion transformation to a copy of the tem- 
plate profile and estimating the induced phase shift be- 
tween the template and its distorted copy. Application of 
the METM-derived instrumental distortion parameters to 
the PSR J1022+1001 template profile yields the system- 
atic timing error shown in Figure [9j The instrumental 
distortion transformation is applied independently in each 
frequency channel before integrating over frequency. The 
induced phase shift is computed using only the total in- 
tensity profile and the Fourier-domai n template matching 
algorithm described by iTavlon (|1992T ) . During the periods 
of maintenance on the Parkes 21-cm Multibeam receiver, 
the solutions derived from observations made with the 
H-OH receiver have been used. Uncalibratcd instrumental 
distortion introduces peak-to-peak variations in system- 
atic timing error of the order of 600 ns; the largest jumps 
in timing error occur when switching between receivers. 
These systematic errors are present in arrival times de- 
rived from either uncalibrated or inaccurately-calibrated 
PSR J1022+1001 observations. As shown in the following 
section, these errors are significantly reduced after calibra- 
tion with the METM-derived solutions. 

3.4. Arrival time estimation and analysis 

The solutions produced using the new METM method 
were applied to c alibrate th e 5-mi nute PSR J1022+1001 
integrations as in lOrd et al.l ([2004). The calibrated data 
were then integrated in frequency and time to form approx- 
imately 64-minute integrations. The observations with 
the greatest S/N were added to form a full-polarization 
template profile with an integration length of ~ 67 hours, 
shown in Figure [TUl 

To enable a quantitative comparison between the anal- 
ysis described in this paper and previous work, the 
PSR J1022+1001 observations were also c alibrated us- 
ing the ideal feed assumption employed by iHotan et al.l 
(20043). Arrival times for the two datasets were then de- 
rived using only the total inte nsity profile ( Tavlor| [l992) 
and matrix template matching (jvan Strated[2006ir These 
variations produced four sets of arrival time estimates, cal- 
ibrated using either full measurement equation template 



matching (METM) or the ideal feed assumption (IFA) and 
timed using either matrix template matching (MTM) or 
standard total intensity (STI) methods. The post-fit ar- 
rival time residuals for these four sets are compared in 
Tabled 

After calibration based on the ideal feed assumption 
(IFA) and arrival time estimation with the standard 
method of template matching using only the total intensity 
profile (STI) the residuals are o f similar quality to those 
presented by IHotan et al.l (|2006f ). That is, compared to 
uncalibrated data, the IFA-calibrated data are equally ad- 
versely affected by distortions to the total intensity profile. 
This is not surprising, given that the IFA do es not apply 
to th e 21-cm Multibeam receiver at Parkes (jvan Stratenl 
2004). Calibration via measurement equation template 
matching (METM) reduces both the standard deviation 
of the arrival time residuals and the reduced \ 2 01 the 
model fit. Assuming that the systematic error removed 
by METM is uncorrelated with the remaining noise in the 
METM-STI residuals, instrumental distortion of the IFA- 
calibrated total intensity profile gives rise to timing errors 
with a standard deviation of approximately 1.3 /is. 

For both METM- and IFA-calibrated data sets, matrix 
template matching (MTM) yields arrival times of similar 
quality. This indicates that MTM is able to correct the 
remaining instrumental calibration errors in the IFA data 
and produce arrival time estimates with reduced system- 
atic error. It is not possible to test MTM on uncalibrated 
data because the variation of instrumental response with 
frequency leads to bandwidth depolarization, which can- 
not be modeled using a single Jones matrix. That is, al- 
though inaccurate, calibration based on the IFA may be 
sufficient to avoid depolarization. 

Focusing on only the METM-calibrated data, compari- 
son of arrival times generated using STI and MTM demon- 
strates that MTM reduces the standard deviation of ar- 
rival time residuals by approximately 35%, which is consis- 
te nt with the re l ative uncertainty predicted in Table [1] and 
by Ivan Stratenl ((2006). When compared to the results of 
the typical data analysis employed in most high-precision 
timing experiments to date (IFA-STI), matrix template 
matching reduces arrival time residuals by over 50%. The 
corresponding improvement in the experimental sensitiv- 
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Fig. 10. — Mean polarization of PSR J1022+1001, plotted as a function of pulse phase using cylindrical coordinates: position angle (top 
panel), linearly polarized intensity, L = \f Q 2 + U 2 (dashed line), circular polarization, V (thin solid line) and total intensity (thick solid 
line). Data were integrated over a 64 MHz band centered at 1341MHz for approximately 67 hours. 



Table 2 

Arrival Time Residual Standard Deviation and Merit 



Method 


c T (M s ) 


X7Wfree 


IFA-STI 


1.9 


2.2 


METM-STI 


1.4 


2.1 


IFA-MTM 


0.92 


1.1 


METM-MTM 


0.88 


1.0 



ity of this dataset is equivalent to that of quadrupling the 
integration length of each observation. 

3.5. Results 

The results presented in this section are based on the 
data calibrated using METM and arrival time estimates 
derived using MTM. Table 13.51 lists the best-fit timing 
model pa rameters obtained using the TEMPQ2 analysis 
software (|Hobbs et alj|2006t lEdwards et al.ll2006t) . These 
results include a significant detection of the Shapiro delay; 
after subtracting two degrees of freedom from the least- 
squares fit with the addition of the shape s — sini and 
range r = parameters, \ 2 IS reduced by ~ 20% (from 
~ 330 to ~ 260). Although the detection is significant, the 
elongated contours of constant A\ 2 shown in Figure [11] 
demonstrate that s and r are highly covariant; the coeffi- 
cient of correlation between these two model parameters 
is —0.99. Furthermore, the constraint on the companion 
mass is not very informative; the 95% confidence interval 
given by the projection of the A\ 2 — 4 contour onto the 
ni2 axis ranges from 0.25 to > 3 solar masses. 

The shape and range parameters are highly covariant 
because, in nearly circular orbits that are only moderately 
inclined with respect to the line of sight, the Shapiro delay 
is readily absorbed in the Roemer delay by modification 
of the Keplerian orbital parameters. The unabsorbed rem- 
nant of the Shapiro delay is best described as a weighted 
sum of harmonics of the orbital frequency using the or- 
thometric parameterization introduced by iFreire fc Wexl 



(|2010[) . Accordingly, Table 13.51 lists the best- fit values of 
the ratio of the amplitudes of successive harmonics 

sini 
S = ; i rr 

1 + | COSl| 

and the amplitude of the third harmonic 
^3 = m 2 S 3 - 

The values of s and r derived from £ and /13 are nearly 
identical to the values yielded by directly modeling the 
shape and range of the Shapiro delay. However, the co- 
efficient of correlation between /13 and <; is only —0.38. 
There is close agreement between the orbital inclina- 
tion angle, i ~ 57°, constrained by the Shapiro de- 
lay and previous estimate s of the inclination of the pul- 
sar s pin axis, C ~ 61° (IXilouris et ail Il998f ) and £ ~ 
55° ( Ra machandran k, Kramerl l2003f), constrained using 
the rotating vector model (Radhak rishnan fc Cookei fl969; 
lEverett fc Weisbe"r3l200l[ ). This agreement is compatible 
with the expectation that the pulsar spin and orbital an- 
gular momentum vectors are aligned during the period of 
mass accret ion that led to the formation of the millisecond 
pulsar fe.g. lThorsett fc Chakrabartvlll999[ ). 

The best-fit model parameters also include the first sig- 
nificant detections of the precession of periastron to and 
the secular variation of the projected semi-major axis of 
the orbit, x. The estimate of Cj ~ (9° ± 3°) x lO^yr" 1 
is consistent with the value predicted by General Relativ- 
ity, Cj ~ 8° x 10 _3 yr _1 , based on the masses of the 
pulsar and its companion yielded by the Shapiro delay 
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Table 3 

Best-fit model parameters for PSR J1022+1001 



Parameter Name 


Value 


Reference clock 


TT(TAI) 


Planetary Ephemeris 


DE414 


P epoch (MJD) 


53869.00016629324189 


Pulse period, P (ms) 


16.4529299518323(2) 


Period derivative, P (1CT 20 ) 


4.33399(6) 


Ecliptic longitude, A (°) 


153.865881(3) 


Ecliptic latitude, /3 (°) 


-0.06(4) 


Proper motion in A (mas yr -1 ) 


-15.97(2) 


Proper motion in f3 (mas yr _1 ) 


38(19) 


Annual parallax, tt (mas) 


1.4(2) 


Binary Model 


DDH 


Orbital period, Pb (days) 


7.805135(1) 


Projected semi-major axis, x (lt-s) 


16.76541(4) 


Orbital eccentricity, e (10~ 5 ) 


9.702(7) 


Epoch of periastron, T (MJD) 


53876.1038(2) 


Longitude of periastron, u> (°) 


97.757(9) 


Orthometric amplitude, /13 (/zs) 


0.52(7) 


Orthometric ratio, <r 


0.5(1) 


Periastron advance, to (°yr _1 ) 


0.009(3) 


Secular variation of x, x (10 -15 ) 


14(1) 



Note. — For each parameter, the formal uncertainty (one standard deviation) in the last digit quoted is given in parentheses. 
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Fig. 11. — Map of A% 2 as a function of binary companion mass and orbital inclination angle. From innermost to outermost, the contours 
represent A\ 2 = 1, 4, and 9. The error bars indicate the best-fit solution and the formal errors (one standard deviation) returned by TEMP02 
for the shape and range of the Shapiro delay. 
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detection. The measured value of x is seven orders of 
magnitude larger than the general relativistic prediction, 
x ~ 5 x 10~ 21 , because it is do minated by the secu- 
lar variation due to proper motion (|Kopeikinlll996l ). This 
geometric effect can be used to place an upper limit on 
the orbital inclination angle; i.e. tani < (ix/x, where 
fi = 42 ± 18 mas yr^ 1 is the total proper motion. The re- 
sulting constraint, i < 83°, rules out only 12% of possible 
orientations and is consistent with the stronger constraint 
yielded by measurement of the Shapiro delay. 

The proper motion of the system also induces a 
quadratic Doppler shift {Shklovskii 1970j) that contributes 
Tfi 2 d/c to observed period derivatives, where T is the pe- 
riod (e.g. spin or orbital), d is the distance to the pulsar, 
and c is the speed of light. The magnitude of the Shklovskii 
effect is estimated using the distance — 750tgo° par- 
sec derived from the annual trigonometric parallax af- 
ter correction for Lu tz-Kelker bias, n = 1.340!q ' jg| mas 
{Verbiest et al.ll2010T ). The Shklovskii contribution to the 
spin period derivative (~ 4.9 x 10~ 20 ) is slightly larger than 
its measured value, which would imply that the intrinsic 
P of the pulsar is negligible. However, the predicted con- 
tribution to the orbital period derivative (~ 2 x 10~ 12 ) is 
around 20 times larger than the upper limit derived by in- 
cluding Pb in the timing model. The relative uncertainties 
in the above estimates of the Shklovskii effect are almost 
100%, primarily due to the large uncertainty in the esti- 
mate of the proper motion in ecliptic latitude. Using only 
the proper motion in ecliptic longitude, the predicted kine- 
matic contribution to Pb is only three times larger than 
the timing-derived upper limit. 

To explain the remaining discrepancy, three other con- 
tributions to the obse rved f\> are considered (as in 
iDamour fc Tavlorlll991[ ). First, the loss of energy due to 
gravitational radiation as predicted by General Relativity 
contributes around —3 x 10~ 16 , approximately 4 orders of 
magnitude smaller than the Shklovskii effect. Second, the 
contribution due to differential rotation about the Galactic 
centre (approximately —7 x 10 -15 ) can also be neglected. 
Finally, the acceleration of the PSR J1022+1001 system in 
the Galactic gravitational potential is considered. Using 
the parallax distance d^ and the Galactic latitude 6 = 51°, 
the component of gravitational acceleration perpe ndicular 
to th e Galactic disk, K z ~ 6.6 x 10~ u m s~ 2 (|Bahcalll 
Il984h contributes approximately —1.2 x 10~ 13 to the ob- 
served P,. At large heights above the Galactic disk, the 
relative uncertainty in K z is more than a factor of two; 
therefore, this contribution is large enough to negate the 
Shklovskii effect, resulting in an observed Pb that is con- 
sistent with zero. 

4. DISCUSSION 

Instrumental distortion of the total intensity profile in- 
troduces significant systematic timing errors that are cor- 
related between pulsars observed with the same instru- 
ment. Therefore, high-fidelity polarimetry is inextricably 
linked with the long-term objectives of high-precision tim- 
ing and the primary goals of pulsar timing array experi- 
ments. This paper presents a novel method of calibrating 
the instrumental response by exploiting the long-term sta- 
bility of the polarized emission from a millisecond pulsar. 
Full-polarization timing observations of PSR J0437— 4715 



are used to model variations in the 21-cm Multibeam and 
H-OH receivers and down-conversion system at Parkes. 
Over a period of ~ 7.2 years, temporal and spectral varia- 
tions of the receptor cross-coupling introduces systematic 
error of the order of 1 [is in arrival time estimates derived 
from observations of PSR J1022+1001. High-fidelity cali- 
bration of instrumental polarization and arrival time esti- 
mation using matrix template matching is demonstrated 
to correct systematic error and double the sensitivity of the 
experiment, yielding significant detections of the Shapiro 
delay, the precession of periastron, and the secular vari- 
ation of the projected semi-major axis of the orbit. The 
improvements in both timing accuracy and precision are 
consistent with the hypothesis that the average polarized 
emission from these millisecond pulsars is stable over the 
timcscales of relevance to this experiment. 

With a modest instrumental bandwidth of 64 MHz and 
median integration length of 64 minutes, PSR J1022+1001 
yields arrival time residuals with an uncertainty-weighted 
standard deviation of only a T = 880 ns, roughly five or- 
ders of magnitude smaller than the pulsar spin period P ~ 
16.5 s. The remarkable timing precision of these ob serva- 
tions supports the conclusions of lHotan et alj (|2004al ). who 
arg ue that the PSR J1022 +100 1 profile variations observed 
by Kramer et all (|l999f ) and iRamachandran fc Kramerl 
(2003|) are not intrinsic to the pulsar and are most likely 
due to instrumental calibration errors. 

The precision of the arrival time data presented in this 
paper places PSR J1022+1001 among the top ten sources 
regularly observed as part of the Parkes Pulsar Timing Ar- 
ray (PPTA) project. Other PPTA sources may also bene- 
fit from the application of METM; e.g. six of the sources 
listed in Table Q] have instrumental distortion susceptibil- 
ity factors ranging from 200 ns to 340 ns. Referring to the 
measured values of the non-orthogonality parameters plot- 
ted in Figure [7l these factors roughly correspond to sys- 
tematic timing error variations of the order of 1 [is. There- 
fore, the timing precision of these pulsars would likely be 
improved by application of the methods developed for this 
study. For PSR J1744-1134 and PSR J2129-5721, ma- 
trix template matching is predicted to yield arrival times 
with greater uncertainty than those derived from the total 
intensity profile. In these two cases, the best results would 
be obtained by calibrating using METM and deriving ar- 
rival times using STL 

In a wider analysis of all of the pulsars in the 
PPTA, Manc hester et al] {2012) calibrated approximately 
5.9 years of PSR J1022+1001 timing observations us- 
ing MEM-derived solutions. Using typical integration 
lengths of 64 min and an instrumental bandwidth of 
256 MHz, this study yielded arrival time residuals with 
an uncertainty- weighted standard deviation of 1.7 fj,s and 
reduced x 2 ~ 9.3. The residual noise is approximately 
2.4 times greater than expected based on simple extrap- 
olation of the METM-STI results presented in Table O 
This discrepancy may be explained by the fact that the 
PDFB instruments have a non-linear response that intro- 
duces over-polarization, which cannot be calibrated using 
MEM. The accuracy of MEM-based calibration is also lim- 
ited by the fact that there is no unique solution to the 
measurement equation when the only constraining trans- 
formation is the geometric rotation of the receptors with 
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respect to the sky (|van Straten! I2004T ). To constrain the 
otherwise degenerate boost along the Stokes V axis, it is 
necessary to include observations of a source that is as- 
sumed to have negligible circular polarization (e.g., for the 
MEM fits performed in ^3.21 observations of Hydra A were 
used to constrain S e ). However, the use of an unpulscd 
source of radiation as a constraint is intrinsically suscep- 
tible to variability in other contributions to the system 
temperature, including receiver noise and ground spillover. 
That is, it is safer to assume that the polarized emission 
from PSR J0437— 4715 is constant than it is to rely on a 
source of unpulsed radiation during c alibra tion. The re- 
sults presented by iManchester et all (|2012[) may also be 
limited by the instability of the receptor cross-coupling, as 
shown in Figure [7J The temporal variations of the non- 
orthogonality parameters are sufficiently resolved only via 
application of the METM method presented in this pa- 
per, which yields a seven-fold increase in the number of 
available calibrator solutions. 

As noted in previous studie s dKramer et all 1999; 
Ra machandran fe Kramer] l2003t iHotan et al.l [2004a; 
Manchester et al1l2012[ ). other phenomena may also cur- 
rently limit the timing precision of PSR J1022+1001. For 
example, owing to its low ecliptic latitude, the radio signal 
from this pulsar is subject to signific a nt dis persive delays 
in the solar wind (|You et al.l l2007bl . 12012ft . No correc- 
tions for dispersion variations were applied to the data 
presented in this work; however, observations made when 
the line of sight to PSR J1022+1001 passes near the Sun 
(around late August of each year) were excluded from the 
data set. In addition to fluctuations in dispersion, the 
observed flux of the pulsar varies as a function of time 
and radio frequency due to scintillation in the interstellar 
medium. The average profile of PSR J1022+1001 varies 
significa ntly as a function of radio frequ ency (e.g. see Fig- 
ure 1 of lRamachandran fc Kram er 2003) and, when modu- 
lated by interstellar scintillation, the frequency-integrated 
mean profile may fluctuate with time. The potentially 
significant arrival time estimation errors induced by this 
effect could be mitigated through the use of a frequency- 
dependent template profile. 

It is reasonable to expect that PSR J1022+1001 tim- 
ing will be improved by the current generation of instru- 
mentation at the Parkes Observatory. The data presented 
in this paper were observed using a system with low dy- 
namic range that performed analog-to-digital conversion 
of the radio signal with only two bits per sample. Two-bit 
quantization is an intrinsically non-linear process and the 
techniques applied to restore linearity (dynamic output 
level setting) and mitigate quantization noise (scattered 
power correction) are based on a linear approximatio n to 
the response of the digitizer (jJenet fc Andersonl lT998). As 
discussed in Appendix [Bj the scattered power correction 
algorithm applied in this analysis is accurate only to first 
order and it is possible that overestimation of the unpolar- 
ized scattered power contributes to the over-polarization 
noted in Figure |5] To overcome the limitations of two- 
bit sampling, a new baseband recording and processing 
system with greater dynamic range was commissioned at 
the Parkes Observatory in 2010 April. Designed in col- 
laboration with the Center for Astronomy Signal Pro- 



cessing and Electronics Research (CASPER) at Berkeley, 
the CASPER-Parkes-Swinburne Recorder (CASPSR) dig- 
itizes a dual-polarization 400 MHz band with eight bits 
per sample and performs real-time radio frequency inter- 
ference excision bas ed on the spectral kurtosis estimator 
(jNita fc Garvl |2010| ). This system currently operates in 
parallel with the third and fourth generations of the Parkes 
Digital Filter Bank (PDFB) instruments as part of the 
PPTA program. 

The calibration techniques applied in this paper could 
potentially be improved by reducing the number of degrees 
of freedom in the estimates of the receptor cross-coupling 
parameters. For example, this could be achieved by fitting 
an analytic model to the temporal and spectral variations 
of the calibrator solutions or by employing a lossy compres- 
sion algorithm that can be applied to irregularly sampled 
data. The smoothing effected by such a transformation 
would reduce the instantaneous noise in the applied cali- 
brator solutions and might also provide a robust means of 
interpolating between solutions. However, further refine- 
ments of the calibration technique may yield only marginal 
improvements to arrival time accuracy and precision. Ta- 
ble[5]demonstrates that, even when the data are calibrated 
using the inaccurate ideal-feed assumption, matrix tem- 
plate matching (MTM) yields arrival time estimates that 
are nearly as good as those derived from data calibrated 
using the new METM technique described in this paper. 
This seems to suggest that, as long as the signal is not de- 
polarized by integrating over time and/or frequency, the 
fidelity of the calibration technique has negligible impact 
on arrival times derived using MTM. 

This conclusion may have a significant impact on the 
design of the next generation of instrumentation for high- 
precision timing. For example, it has been asserted that 
in order to achieve timing precision of the order of 100 ns, 
the Square Kilometre Array (SKA) must ac hieve net po- 
lariza tion purity corresponding to /3 < 10~ 4 (ICordes et al.l 
|2004|) . However, even when non-orthogonality as large as 
(3 ~ 10 -2 is ignored, as is the case when the 21-cm Multi- 
beam receiver at Parkes is assumed to be ideal, MTM 
yields accurate arrival time estimates. That is, by relaxing 
the requirement for polarization purity, the application of 
MTM has the potential to reduce the cost of SKA devel- 
opment for high-precision timing. 

All of the software required to perform measurement 
equation modeling, polarimetric calibration, and ma- 
trix templ ate matching is fr eely available as part of 
psrchive (IHotan et alj|2004bl): the use of this software is 
demonstrated by Ivan Straten et all (|2012h and more fully 
documented onlin^t 
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Mike Keith, Michael Kramer, Stefan Oslowski, and John 
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APPENDIX 

PROBABILITY DISTRIBUTION OF THE POLARIMETRIC INVARIANT 

The invariant profile is formed by computing the Lorentz interval S as a function of pulse longitude, where S 2 = S 2 — \S\ 2 . 
Assuming that the noise in each of the Stokes parameters is independent and normally distributed with standard deviation 
?, define the normalized Stokes parameters S' k = S^-A- The noise in Sq has a noncentral x 2 distribution with one degree 
of freedom and noncentrality parameter A = S^ 2 . Likewise, the noise in \S'\ 2 has a noncentral x 2 distribution with three 
degrees of freedom and A = IS 1 '! 2 . The distribution of the noise in S 2 is given by the cross-correlation of the distributions 
of S'q and IS 1 '! 2 ; examples are shown in Figure [A121 The variance of S 2 is equal to 4||S"|| 2 + 8, where \\S\\ is the 
Euclidean norm of the Stokes four-vector, defined by HS^ 2 = S 2 + \S\ 2 . As all four Stokes parameters vary as a function 
of pulse longitude, the noise in S" 2 is heteroscedastic, which violates one of the basic premises of most template-matching 
algorithms. 

Both the distribution of the invariant interval and its tendency toward zero for highly polarized radiation adversely 
affect its usefulness in template matchi ng (as a replaceme nt for the total intensity) and as a normalization factor during 
measurement equation modeling (MEM; Ivan Strat cn 2004). Gain normalization is necessitated by interstellar scintillation, 
which causes the received flux of the pulsar to vary as a function of time and radio frequency. To address the problems 
with the invariant interval during the application of MEM in i j3.2[ the Stokes parameters were normalized by the sum 
of the invariant interval over all on-pulse longitudes. For the typical observation of PSR J0437— 4715, this summation 
reduces the relative error of the normalization factor by ~ 97%, thereby yielding more normally distributed normalized 
Stokes parameters. 

FOLDED PROFILE SCATTERED POWER CORRECTION 

The data presented in this paper were processed using a new scattered power correction algorithm that can be used to 
correct digitization distortions to folded pulsar profiles. In the case of two-bit sampled voltages, the scattered power (or 
quantization noise) is ~ 1 2% of the total power. T he derivation of the method starts with the mean digitized power a 2 
given by equation (A5) of lJenet fc Anderson! (jl998l hereafter JA98), 

a 2 = 5>($)/($), (Bl) 

where <& is the fraction of samples that fall between the chosen thresholds, /($) is the digitized power as a function 
of <£>, given by equation (A4) of JA98, and V{&) is the discrete probability distribution for $, given by equation (A6) 
of JA98. The parameter $ is eliminated by the summation in the above equation; however, both /($) and V^) are 
also parameterized by the expectation value (<£>). Therefore, equation (|B1[) represents the relationship between the mean 
digitized power and the mean value of $. That is, given the mean digitized power er 2 , equation (|B1[) can be inverted to 
compute ($), which can in turn be used to estimate the mean undigitized power a 2 and the mean scattered power 1 — A 
via equations (45) and (43) of JA98, respectively. 

Equation (A5) of JA98 can be inverted using the Newton- Raphson method and the partial derivatives of equations (A4) 
through (A6) of JA98 with respect to (<&). These are simplified in the case of two-bit sampling by noting that equation 
(A4) of JA98 reduces to 

/($) = ($)y ; 2 ($) + (l-($))y 2 ($) 
where 2/3(3?) and ua(^) are given by equations (41) and (40) of JA98. 




Fig. A12. — Probability density of the centralized polarimetric invariant, S 2 = I 2 — P 2 — A, where / is the total intensity, P is the polarized 
flux, and A = {I) 2 — (P) 2 is the noncentrality of S . The standard deviation in each of the four Stokes parameters is unity; therefore, the 
mean value of S is —2. The four curves correspond to {/} = (solid), 1, (long dash), 2 (short dashed), and 3 (dotted). In each case, (P) = 
0. For large values of {/), the probability density of S 2 approaches a normal distribution. 
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The folded profile scattered power correction algorithm is based on the following assumptions and approximations. 
First, at the time of recording the astronomical signal, it is necessary that the baseband voltages are sampled using the 
optimum two-bit input thresholds defined in Table 1 of JA98. To first order, this condition is satisfied by excluding data 
with excessive two-bit distortion as defined by equation ([9]) of this paper. Second, equations (45) and (43) of JA98 only 
approximately relate the expectation values of <£>, cr 2 and A. For example, the relationship between A and a 1 defined by 
equation (43) of JA98 is concave down (see Figure 3 of JA98); therefore, Jensen's inequality dictates that the value of 
A estimated from the mean undigitized power will be greater than the expectation value of A. Consequently, the mean 
fractional scattered power may be systematically underestimated. This limitation cannot be overcome without prior 
knowledge of the distribution of total intensity fluctuations intrinsic to the pulsar signal. Finally, it is assumed that the 
signal in each frequency channel has not been significantly altered, such that there is a well-defined relationship between 
the mean undigitized power and the mean digitized power. For example, after phase-coherent dedispersion removal, the 
flux density over a given time interval no longer represents the voltage fluctuations in the digitizer over that interval; the 
previously smeared pulsar signal will be recovered and the scattered power will be smeared. To estimate the scattered 
power using coherently-dedispersed digitized power, the dispersion smearing in each frequency channel must be less than 
or of the order of the time resolution of the folded profile. This condition is satisfied in the 20-cm timing observations of 
PSR J0437-4715 and PSR J1022+1001 presented in this paper. 
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